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Abstract 



Two methods are compared that are used in path integral simulations. Both 
methods aim to achieve faster convergence to the quantum limit than the so- 
called primitive algorithm (PA). One method, originally proposed by Taka- 
hashi and Imada, is based on a higher-order approximation (HOA) of the 
' quantum mechanical density operator. The other method is based upon an 

^ ■ effective propagator (EPr). This propagator is constructed such that it pro- 

QQ , duces correctly one and two-particle imaginary time correlation functions in 

', the limit of small densities even for finite Trotter numbers P. We discuss the 

conceptual differences between both methods and compare the convergence 
rate of both approaches. While the HOA method converges faster than the 
EPr approach, EPr gives surprisingly good estimates of thermal quantities 
already for P = 1. Despite a significant improvement with respect to PA, 
\ neither HOA nor EPr overcomes the need to increase P linearly with inverse 

Q I temperature. We also derive the proper estimator for radial distribution func- 

^ ' tions for HOA based path integral simulations. 



I. INTRODUCTION 

Path integral Monte Carlo (PIMC) [|l| and path integral molecular dynamics (PIMD) 
have proven useful in the atomistic simulation of quantum effects occuring in condensed mat- 
ter at low temperatures. The broad range of applications includes among other things super- 
fluid ''He [0], isotope effects in crystalline rare gas solids as well as phase transitions with 
strong quantum mechanical effects in Josephson junctions and molecular solids 0,0. 
Recently, path integral methods have also been applied to calculate low-temperature mate- 
rial properties of systems that are computationally less easily tractable than Lennard- Jones 
type systems: solid germanium [|10|, crystalline polyethylene |jTl|, diamond 0], silica |]13 



and even Wigner crystals [|I^, to name a few. 

A disadvantage of PIMC and PIMD is the increase of necessary computing time tcpu with 
decreasing temperature T at a given required accuracy. Using the so-called primitive ap- 
proximation together with the most efficient sampling algorithms that completely eliminate 
critical slowing down (see Ref. P for a thorough discussion), it is not possible to overcome 
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^cpu oc 1/T. Therefore different improvements on the primitive algorithms have been sug- 
gested that render path- integral simulations more efficient. One may subdivide the bulk 
of such attempts into three categories (i) methods that are based on higher-order approx- 
imants (HOA) of the high-temperature density matrix |15|, (ii) methods that use effective 



propagators (EPr) that automatically yield the proper two-particle behavior PJT6|] , and (iii) 
methods that decompose the Hamiltonian into a harmonic and an anharmonic contribution 
before applying the Trotter formula |T^,|TB[. The latter category, also refered to as the effec- 



tive potential (EP) method, will not be discussed any further in this study, mainly due to its 
unfavorable scaling of tcpu with particle number N. Additionally, harmonic approximations 
are problematic for many systems of interest, in particular those involving '^He and ^He. Yet, 
another advantage of HOA and EPr methods over the EP method is that the pathological 
behaviour of the attractive Coulomb potential is overcome automatically |T6| , p!9| , |20[] . 

In this study, we want to compare the convergence of the HOA method, the EPr method, 
and the primitive algorithm for a simple model system. As long as the interaction potentials 
are well-behaved, the convergence does generally not depend on the specific form of the 
potential ||2T[. The test model system should of course be chosen such that it does not favor 
intrinsically one approach over the other, e.g., we may not chose a two-particle system, 
because then the EPr model would be exact per definition. As we are not interested in 
the EP approach, we can chose a simple monoatomic chain with harmonic next-neighbor 
coupling. This choice of system enables us to do the bulk of the calculations analytically so 
that statistical error bars are eliminated completely. We also want to investigate how the 
Trotter number P necessary to keep the systematic errors below a well-defined percentage 
scales with inverse temperature for the different approaches investigated in this study. 

It should be emphasized that the HOA method and the EPr method are conceptually 
different. In an EPr path integral simulation, one tries to generate radial distribution func- 
tions that are in the quantum limit (at least in a low-density approximation). Evaluating 
observables such as the thermal expectation value of the potential energy (Ipot) is done by 
simply evaluating the operator of the potential energy at the given distance. In a HOA path 
integral simulation, generalized estimators have to be defined even for those observables that 
are orthogonal in real space. We will comment on this in more depth in the following section. 
In particular we will derive an expression for the HOA estimator of the radial distribution 
function, which has not been given hitherto, which might explain the sparse use of the 
method in the literature. The new estimator will be used to calculate the argon argon radial 
distribution function of those atoms in a simple three-dimensional Lennard- Jones crystal. 
The different methods employed in this study will be outlined in Sec. II. In this section we 
will also present a simplified approach to the EPr approach, which we call reduced effective 
propagator (r-EPr) approach. The results will be presented in Sec. Ill, and a short summary 
is given in Sec. IV. 

II. METHODS 

A. Primitive Algorithm 

The primitive algorithm for path integrals is based on Feynman's idea to represent the 
partition function of a quantum mechanical point particle Z{/3) as a partition function of a 
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classical ring polymer p2| , p3| . The potential V^p({?"}) of the classical ring polymer has the 
form: 



Vr,{{r}) = E 



t=i 



2in-rt+i) +V{rt) 



(1) 



with (3 = l/Zc^T. Tt represents the position of monomer t in the ring polymer (r^ = r^+p), 
and V is the real (physical) potential, t is sometimes interpreted as an imaginary time and P 
is commonly called the Trotter number. The PIMC or PIMD program are then assumed to 
generate distributions such that the probability of configuration {r} to occur is proportional 
to exp[— /5Vrp({r})/P]. All thermal expectation values of observables diagonal in real space 
can be determined directly from the configurations, e.g., 



M P 



{V) 



lim lim 

P^oo M^oc 



MP 



(2) 



=1 t=l 



where rt^i is the position of the t'th monomer in the i'th Monte Carlo step and M is the 
number of observations in the Monte Carlo simulation. We refer to Refs. for further 
details on the primitive algorithm. 



B. Higher-Order Approximant Method 



The HOA method is based on a fourth-order Hermitian Trotter decomposition of the 
high-temperature density matrix p3. The decomposition was first applied to continuous 



degrees of freedom by Takahashi and Imada |jT5[ as well as by Li and Broughton |T9[ . The 
basic idea of the decomposition is to approximate the high-temperature density matrix 
p = exp{-f3H/P) with 



P 



^ g-/3V'/2Pg-/3T/2Pg-/3V'cor/Pg-/3T/2Pg-/3y/2P 



(3) 



where H = T + V corresponds to the Hamiltonian and Kor = /3^[[V", T], V"]/24P^ is a 
correction term. T and V are usually chosen to be kinetic and potential energy, respectively. 
For this decomposition, the correction energy V^or can be written as 



/o2fc2 N 1 



(4) 



where m„ corresponds to the mass of particle n. Owing to the temperature-dependence of 
Vcor, thermal expectation values of observables O have to be reevaluated with respect to the 
primitive algorithm, e.g., averages of functions diagonal in real space read 



1 ^ i 

TTsli^ EE 0(W) 



MP 

N 

+ E 
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(5) 
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where V on the right-hand side of Eq. ^ only includes the original potential and not the cor- 
rection term V^or- The accuracy of the HOA method outlined above allows one to determine 
thermal expectation values with a leading correction of while the primitive algorithm 

has leading corrections in the order of [[TStpiQipI | . 



Eq. (^) allows one easily to find the estimator for the potential energy to be \^ + 2V^or, 
see Refs. [p!5l , |l9| for further details on the calculation of thermal expectation values. To 



our knowledge, however, it is has not yet been discussed that even the estimator for radial 
distribution functions g{r) needs to be altered with respect to the primitive approach, for 
which the estimator can be written as 

g;':^ir)^5ir-\rt,i,n-rt,,n' \) / r\ (6) 

rt,i,n denoting the position of particle n of the fs monomer in ring polymer (particle) n. 
Applying Eq. (^) to Eq. (^) leads to a shift of the estimator for the distance between particle 
n and n'. Simply applying Eq. (|^) to the operator for the square of the distance between 
particle n and n' yields the estimator r^*^^™ for the distance between particle n and n', which 
is found to be 

Thus r^n^ should replace r in the argument of the 5-function on the right hand side of 
Eq. (^) in order to calculate g{r). Based on this relation, one may say that the HOA 
method is not an importance sampling algorithm in the sense that the probability for a 
configuration to occur in the simulation is proportional to the diagonal elements of p in a 
real space representation. 



C. Effective Propagator Approach 

An alternative non-primitive method to the HOA method is to construct effective po- 
tentials such that the two-particle propagators are reflected accurately, e.g., it is correct in 
all orders of h. This effective potential is then used in a multi-particle simulation and hence 
produces the proper thermal behaviour in the low-density limit. No effective estimators 
have to be defined for observables diagonal in real space and in this sense, the approach is 
an importance sampling algorithm. The effective propagator (EPr) method is particularly 
useful for ill-behaved potentials such as the attractive Coulomb potential |T^. One may 
write the high-temperature two-particle density operator in the following way: 

(rir2 I e~^^/^ \ r'ir'2) oc exp 



(8) 

Kfj is a function that depends on ri,r'^,r2, and r'2. Therefore the interaction can be said 
to be non-local in imaginary time, e.g., in the primitive decomposition ri does not couple 



Ik, 



X exp 



(3 mP^ 



{(ri-r;)^ 



+ r2 
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directly to r'g. The calculation of both the diagonal and the non-diagonal elements in more 
than one dimension is not trivial for non-harmonic potentials and we refer to Ref. for an 
in-depth discussion of that problem. For our one-dimensional model system, however, the 
approach can be simplified significantly, i.e., it can be solved analytically up to a summation 
over finite number of terms. This will be done in the Section fT\\ . 



1. The r-EPr method 

In general the implementation of the full two-particle pair propagator with correct di- 
agonal and non-diagonal elements is difficult because, one has to use the two-particle high- 
temperature density matrix (HTDM) (rir2 | e~^^^^ \ r'ir'2). In principle one can evaluate 
this expression prior to the simulation. However, this is a difficult task due to the dimen- 
sionality of the matrix, in particular for anharmonic d = 3 dimensional systems. Even if 
one expresses the HTDM in the center of mass system of the coordinates ri, r2, r'l, r'2 and 
even if one uses an appropriate orientation of the axes within a molecular-fixed frame, one 
is a left with a three dimensional table and transformations between molecular fixed frames 
and laboratory system. It is also difficult to find good fit functions for the effective potential 

The main idea of a reduced effective potential approach (r-EPr) is to only incorporate 
those corrections to the primitive decomposition that are local, e.g., those corrections that 
involve coordinates at different Trotter indices are neglected. In other words, if one carries 
out a simulation at inverse temperature (3 with Trotter number P, one uses an effectively 
classical potential that reproduces the correct two-particle distribution function in the low- 
density limit at inverse temperature (3/P. This approach is similar to an idea suggested 
by Thirumalai et al. |2^, who constructed effective interaction potentials from the diagonal 



elements of the high-temperature density matrix, see also a related paper by Pollock and 
Ceperley [p6 . 



III. RESULTS 

A. Linear, harmonic chain 

In order to analyze the convergence of primitive algorithm, HOA and EPr method re- 
spectively, we chose a one-dimensional linear chain with harmonic next neighbor coupling 

1 ^ 

V = -Y.k{rn-rn^ir. (9) 

^ n=l 

Periodic boundary conditions, r^+i = ti are applied and the masses m are identical for all 
atoms. 

HOA and EPr invoke correction terms in the potential energy of the ring polymers with 
respect to the original expression of the primitive approach that is given in Eq. (|I|) for a one 
particle problem. All three approaches can be represented as the limiting case of a (1+1)- 
dimensional solid with harmonic coupling between nearest neighbors, next nearest and next 
next nearest neighbors. A graphical illustration is given in Figure |I|. 
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The new effective energy V^p that enters the Boltzmann factor reads: 

N P 



^ n=l t=l ± 

+ ~k2{rt,n - n,n+2f ] . (lO) 



The expression can be diagonahzed in terms of Fourier modes 

p N 

EE ^t.ne^^"'"^^^^'^"*^^- (11) 



'rp 

r~\ p N 



t=l n=l 



It is convenient to introduce k' = k + k\ and k' = n + ki with n = mP'^ /{[3'^% ) in order to 
reexpress V^p as 



rp 

P N 

2 



^rp ~ O E E ^'J-'^ I 1^ (12) 



with 



aj=l (jr=l 



fcg,<^ = K'4sin^ (^^) ^'^si'^^ 



[sin^ (vr(g/A^ + cu/P) + sin^ {TT{q/N - uj / P) 
+Msin2(^^2g) . (13) 

The partition function Zc for the classical system illustrated in Fig. |1] can then be reduced 
to NP Gaussian integrals. Z^. is proportional to 

N P 



oc n n ■ (14) 

q=l LJ=1 

For the different algorithms we find different functions for ki,k2,iii, and K2. The expressions 
for these effective coupling coefficients are summarized in Table |. The expressions A and 
C used in Table | for the EPr treatment will be given and derived below. 

We will now calculate the thermal expectation value of the potential energy (Vp) for a 
given Trotter number P. One can expect that the convergence rate does not depend on the 
observable, which is why it is sufficient to only investigate (Vp). The calculations for the 
three approaches will be separated into three subchapters. 

1. Solution for the primitive Method 

In the primitive method, the usual thermodynamic relationships can be used in a straight- 
forward way without modification, albeit their use can be impractical for technical reasons. 
In order to calculate the thermal expectation value of (Vp), one can use the relation: 
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where all effective coupling coefficients vanish. This relationship simply follows from the 
formal expression for the quantum mechanical partition function of a linear, monoatomic 
harmonic chain. The final expression for (Vp) then is: 

(M^^EEi^H^. (16) 

This equation can be interpreted as follows: Equipartition requires that kq ^^j | Xq^^j p — kj^T. 
Hence kq^^^ appears in the denominator. The amount of real potential energy in these modes, 
however, is only kq \ Xq^^^ P /2, where kq is the stiffnes that can be associated with the mode 
g in a classical linear chain. 



2. Solution for the HOA method 

Eq. (^) also holds for the HOA method, because the only modification with respect 
to PA is that a better approximant for the high-temperature density matrix is emploied. 
However, the correct coupling coefficients ki and k2 have to be used. The expressions for 
ki and /c2 stated in Table |I| are obtained in a straightforward way by inserting Eq. into 
Eq. (H). This leads to the expression 



k ^ P { ( 



■~f4V*-=f^2^)IA,T. (17) 



3m \PJ \N 

where k^^^ refers to that expression for fcg in Eq. (|T^) which is obtained by inserting the 
HOA values for ki and k2. The same result for (Vp) could have been obtained by calculating 
the second moments of the eigenmodes (| Xq^^ p) from equipartition. The resulting (| Xq^^ p) 
could then have been used to calculate the proper HOA potential energy estimator V + 2Vcot- 



3. Solution for the EPr method 

In the EPr approach, does not follow from a Hermitian decomposition of the high- 
temperature density matrix. Therefore Eq. (p3D can not be emploied for the calculation of 
(Vp). However, we can use the fact that EPr is an importance sampling method. Since we 
can obtain (| Xg ^^ p) from equipartition and since the real potential energy of mode g in a 
(classical) linear chain is proportional to A; sin^(7rg/A^), it is possible to use Eq. (|T^ where 
the denominator kq^^) is taken from the EPr column in Table |. Hence, 

q=lw=l l^q^u 
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We are now concerned with the derivation of the expressions for A and C in Table |. 
To do this, one has to consider a dimer in which the two atoms are coupled by a harmonic 
spring of stiffness k. One then transforms the dimer Hamiltonian into the center-of-mass 
system with R = {ri + T2)/2 the center-of-mass coordinate and AR = (ri — the reduced 
distance between the two atoms. These two coordinates can be associated with two modes, 
the center of mass mode with mass M = 2m and the internal oscillator mode with stiffnes k 
and reduced mass fi = m/2. Both the free particle density matrix and the high-temperatur 
density matrix (HTDM)of an oscillator with spring constant k and mass /i are known exactly. 
The free particle HTDM is simply proportional to exp[— MP(i? — R'Y /2l3fi?] (see Eq. (|I|)) 
while the oscillator's HTDM is given by E 



p(AR, AR',p/P) 

j —y/JIk 

''''^[2hsmh{f) 



X 



\\ 2'Khsinh{2f) 
[AR^ + AR'^) cosh(/) - 2ARAR' 



(19) 



The prefactor on the right-hand side of Eq. ([I9[) provides an irrelevant offset in V^s, which 
will be neglected in the following treatment. One then needs to transform the product of the 
free HTDM and the internal oscillator HTDM back into the laboratory system and express 
the effective potential according to Eq. (H). With the definition of 



we obtain the parameters A and C 



/2 1 sinh(/) 



(20) 

(21) 
(22) 



that were introduced in Table |. 



4- Solution for the r-EPr approach 

According to Sec. [II C 1| , we need to know the quantum-mechanical radial distribution 
function g{r) at inverse temperature /5/P for the reduced harmonic oscillator in order to 
construct the r-EPr effective potential V^g. For this purpose, it is sufficient to know the 
internal energy of the reduced harmonic oscillator, since g{r) is a simple Gaussian at all 
temperatures. Hence we need to find the effective coupling constant k' which generates a 
second moment (x^) in a classical treatment at temperature PT such that the the real poten- 
tial energy k{x'^)/2 corresponds to the quantum limit at temperature PT. This condition, 
which defines A;', can be written as 

{V{NP))^^^^, = \^kBTP, (23) 
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where (^(/3/P))exact is the correct thermal potential energy for an oscillator at temperature 
PT. For a harmonmic oscillator, {V {(3 / P)) exact is half the internal energy U, given by 
U = 0.5 hu coth(/5/ico'/2) with frequency u = \Jk/ ft = ^2k/m. We solve for k' and find 

k' = Ak. The parameter A is given in Eq. (^) and ki = k{A — 1) as stated in Table | 
follows. 



5. Comparison of the methods 

The main issue of this study is the analysis of the convergence of thermal expectation 
values such as the potential energy (Vp) to the proper quantum limit as a function of Trotter 
number P. We consider a linear chain consisting of = 5 atoms and periodic boundary 
conditions. The convergence does not depend on in a qualitative way. It is examined 
at a fixed thermal energy well below the Debye frequency of the chain, namely at inverse 
temperature /? = 64/ {h^ k/m). A linear plot of (Vp) is shown in Fig. |^. 

It can be seen that at P = 1 the EPr and the r-EPr method start off with estimates that 
are very close to the quantum limit while PA and HOA start off with an estimate near the 
classical value. Upon increasing P the EPr approaches the proper value from below, while 
for the r-EPr method the deviation between estimate and proper result first increases before 
it decreases again. At a Trotter number P 64, the HOA method becomes similarly good 
as the EPr approach. In order to address the convergence in a more quantitative way, it is 
convenient to analyze the relative deviation of (Vp) from the exact value (Vg) as a function 
of P in a double logarithmic plot, which is done in Fig. |^. 

It can be seen that at large Trotter numbers the HOA method converges with P~^ to 
the quantum limit, while all other methods only converge with P~^. The prefactor for the 
EPr, however, is much smaller than for r-EPr and PA. The value of P at which convergence 
starts is similar in all approaches, e.g., ksTP is larger than but in the order of h^k/m. The 
accuracy where HOA becomes better than EPr is 1.7%. For this particular model system, 
this value of 1.7% was found to be independent of temperature. We expect it to be similar for 
all systems that are dominated by harmonic interactions. We will now address the question 
how we have to increase P for the various approaches if we lower T and require the relative 
accuracy to be constant, e.g., 1%. The results are shown in Fig. H. 

In order for the relative error to be constant, all methods require that P increases linearly 
with inverse temperature (5. The r-EPr approach, which is a little more difficult to implement 
than the PA method, requires slightly reduced Trotter numbers with respect to PA. We 
want to emphasize that the behavior shown in Fig. ^ is qualitatively similar if the accuracy 
criterion for P is changed, however, the stricter the criterion the larger the gap between EPr 
and HOA. This trend can be seen in Fig. ^), where we require one time 0.1 % accuracy 
instead of 1 % as shown in Fig. ^a). Only if one is confined to the use of very small P, EPr 
might be the better choice. One may conclude that the optimal method depends on the 
desired accuracy. 
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B. Crystalline Argon 



It has been pointed out in Section [11 B| that in order to calculate radial correlation 
functions, it is necessary to correct the estimator for g{r). An important question to address 
is how well the HOA approach allows one to calculate g{r). In order to examine this issue 
we apply the HOA method (making use of the proper estimator, see Eq. (0)) to crystalline 
argon. The results are presented in Fig. ^ 

It is interesting to note that g{r) is too broadened for the HOA approach while it is too 
narrow using the PA algorithm. Obviously, the agreement of the P = 12 HOA calculation is 
already very close to the quantum limit. This is a little surprising as the product of ksTP is 
still far below the thermal energy of the Debye temperature, which is about To 70 K. We 
want to note that if the g{r) are obtained without corrections, the agreement is distinctly 
reduced. This might be the reason why the convergence of g{r) reported by Li and Broughton 
for the attractive Coulomb potential was so slow In our case, omitting the corrections 



to the g{r) estimator leads to peaks in g{r) that are even sharper than those obtained with 
PA using identical Trotter numbers. 



IV. SUMMARY 

In this study we have compared the convergence to the quantum limit for different path 
integral approaches. As the convergence rate does not depend on the specific model system 
(as long as the potentials are well-behaved), we have focused our attention to a linear 
chain of harmonically coupled atoms. We disregarded methods making explicit use of a 
decomposition into a harmonic and an anharmonic part of the Hamiltonian among other 
reasons due to the unfavorable scaling of the numerical effort with system size. 

The three main approaches investigated were the so-called primitive algorithm (PA), a 
method based on a higher-order approximant (HOA), and an approach in which an effec- 
tive potential is constructed such that the one and two-particle high-temperature density 
matrices are reproduced exactly in the limit of small densities. We called this approach 
the effective propagator (EPr) method. The latter approach can be further reduced by only 
taking into account corrections that are local in imaginary time leading to the r-EPr method. 

We emphasized the different spirit of EPr and HOA methods: In the EPr method, 
observables orthogonal in real space can be evaluated directly, whereas the HOA method 
requires effective estimators. In this sense, HOA can be seen as a non-importance sampling 
technique that requires the use of effective estimators even for radial distribution functions. 
The need to alter estimators for radial distribution functions obtained in simulations taking 
into account quantum effects in an effective potential had not yet been discussed hitherto. 
We showed that the proper estimators lead to a very fast convergence to the quantum limit 
in a simple argon crystal. 

For the linear monoatomic chain investigated in this study, we found that the correc- 
tions based on HOA vanish with while all other methods, PA, EPr, and r-EPr, have 
corrections in the order of P^^. The prefactor is similar for PA and r-EPr and distinctly 
smaller for EPr. Both EPr and r-EPr, however, give rather accurate results at small P, e.g., 
the relative error in the thermal expectation value of the potential energy is smaller than 
10% in those approaches, while PA and HOA differ by nearly 100% at low temperatures and 
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P — 1. The Trotter number where convergence starts is similar in all approaches. If one 
requires high accuracy, e.g. 1%, one needs to increase the Trotter number P linearly in all 
approaches. 

What can we conclude for path integral simulations? The use of HOA is certainly a 
little more (about twice) CPU time expensive than that of EPr, r-EPr, and PA (which all 
require approximataly the same amount of computing). This is because in an HOA based 
simulation, we need derivatives of the interaction potential that are one order higher than 
those in simulations based on EPr, r-EPr, and PA. However, one is rewarded with the best 
convergence to the quantum limit in a HOA based path integral simulation. While EPr 
also results in a significant improvement with respect to PA, it is plagued with a tedious, 
non-trivial implementation procedure. Thus HOA should be the method of choice unless 
many-body correlations do not play a significant role like in gaseous helium. 
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TABLES 





PA 


HOA 


EPr 


r-EPr 
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k(A - 1) 
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{hl3kf/{l2mP^) 
















-kC 





^2 








\kC 






TABLE L Expressions for the effective coupling coefficients that are 



represented in Fig. || 
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FIGURES 
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FIG. 1. Illustration of coupling between atoms. Atoms can only os- 
cillate in horizontal direction. Vertical springs also act in horizontal 
direction. The straight horizontal lines represent springs between near- 
est neighbors of stiffness k, the solid vertical lines springs of stiffness 
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FIG. 2. Thermal expectation value of the potential energy (Vp) as 



a function of Trotter number P; j3h\Jk/m = 64. 
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FIG. 3. Relative error of the potential energy for an A?^ = 5 chain at 
PUu = 64 as a function of Trotter number P. 
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FIG. 5. Pair correlation function g{r) of crystalline argon at T = 2K 
calculated with the PA and HOA algorithm for Trotter numbers P = 8 
and P = 12. As a reference a quasi-exact correlation function (HOA, 
P = 256) is included. 
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